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ABSTRACT 

Plane viscous channel flows are perturbed and the ensuing initial-value 
problems are investigated in detail. Unlike traditional methods where trav- 
eling wave normal modes are assumed for solution, this work offers a means 
whereby completely arbitrary initial input can be specified without having 
to resort to eigenfunction expansions. The full temporal behavior, including 
both early time transients and the long time asymptotics, can be determined 
for any initial disturbance. Effects of three-dimensionality can be assessed. 
The bases for the analysis are: (a) linearization of the governing equations; 
(b) Fourier decomposition in the spanwise and streamwise directions of the 
flow and; (c) direct numerical integration of the resulting partial differen- 
tial equations. All of the stability data that are known for such flows can 
be reproduced. Also, the optimal initial conditions can be determined in 
a straight forward manner and such optimal conditions clearly reflect tran- 
sient growth data that is easily determined by a rational choice of a basis for 
the initial conditions. Although there can be significant transient growth for 
subcritical values of the Reynolds number using this approach, it does not 
appear possible that arbitrary initial conditions will lead to the exception- 
ally large transient amplitudes that have been determined by optimization 
of normal modes. The approach is general and can be applied to other 
classes of problems where only a finite discrete spectrum exists, such as the 
boundary layer for example. 
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for Computer Applications in Science and Engineering (ICASE), NASA Langley Research, 
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1 INTRODUCTION 


The subject dealing with the fate of disturbances in well established 
flows remains a major topic of fluid mechanics. As the means of obtaining 
quantitative results for any specific problem has improved, the number of 
contributions in this area has grown almost without bounds. And, indeed, 
almost every conceivable prototypical flow has been subjected to scrutiny 
in this manner. The scheme is well conceived: introduce perturbations into 
a defined basic flow, linearize the governing equations, and then determine 
from the initial-value problem the resulting dynamics. In principle this can 
be done but, in practice, it is a formidable task. Some of the difficulties entail 
non self-adjoint and singular differential equations among other obstacles. 
Even the use of the computer and numerical methods have not been easily 
adaptable. Thus, the net result is that (almost without exception) it was 
considered adequate if a flow was determined to be stable or unstable. This 
was done by assuming a separable normal mode solution in the form of 
traveling waves and then establishing the existence of at least one unstable 
eigenvalue. No attention was directed to any particular initial- value or the 
transient period of the dynamics. Indeed, this part of the evolution was not 
thought to have any significance and, in view of the complications involved 
in the linear ordinary differential equations, it was left to speculation. 

More recently the early transient period for the perturbations has been 
shown to reveal that a superposition of decaying normal modes may grow 
initially albeit decay as time goes on. Although the basic origins and recog- 
nition of this type of response should properly be given to Kelvin (1887) and 
Orr (1907a, 1907b), there is an ever widening probe being made; cf. Boberg 
and Brosa (1988), Gustavsson (1991), Butler and Farrell (1992), Reddy and 
Henningson (1993), or Trefethen et al. (1993) as major references. Briefly, it 
has been shown that transient growth can be significant even for subcritical 
values of the Reynolds number and therefore nonlinearity may ensue, mak- 
ing an exponentially growing normal mode a moot point. In point of fact, 
this sort of behavior is not hard to grasp because any non self-adjoint differ- 
ential equation will have eigenvectors that are not orthogonal and therefore 
there can be algebraic behavior for early time. It is clear that the least 
damped normal modes should be the ones that are dominant during this 
early period. 

Other approaches for demonstrating algebraic growth have been given 
by Gustavsson (1979), Benney and Gustavsson (1981), and Criminale and 
Drazin (1990). In the first case, Laplace transforms were used for the distur- 
bance equations and algebraic behavior was found because there are branch 
cuts needed for the inversion of the transforms; poles correspond to expo- 
nential growth or decay. In the second work, it was shown that, if three- 
dimensional disturbances are considered, then there can be a resonance be- 
tween the normal modes of the Orr-Sommerfeld equation and those of the 
Squire equation. Such resonance does occur and it is for damped expo- 
nential modes. The last paper dealt with the existence of the continuous 
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spectrum as well as the discrete normal modes but, again, it is algebraic 
temporal behavior that results. The Laplace inversions of Gustavsson that 
led to algebraic growth were also due to the continuous rather than the dis- 
crete spectrum but, unlike the Criminale and Drazin presentation, general 
initial- values were not taken into account. 

The results of the various investigations that have been made for channel 
flows suggest that the early transient growth can be extremely large (factors 
of 20 to 1000) at subcritical values of the Reynolds number when consid- 
ering optimal growth determined from variational techniques. It should be 
emphasized that not one exponentially growing mode is involved in this 
process and, although not explicitly stated, the implication is that no mat- 
ter what may be the initial input, this exceptional growth will be present. 
Such a premise is subject to scrutiny. Granted variational calculations will 
assign values to the coefficients of the modes of the expansion in terms of 
the eigenfunctions, but such a determination may not be the same as for a 
specific initial- value distribution used to determine the relative values of the 
coefficients of the eigenmodes. As outlined in Betchov and Criminale (1967) 
and Drazin and Reid (1981), the use of the eigenfunctions in the initial- value 
problem requires the adjoint differential equation in order to fix the coeffi- 
cients. Moreover, it is known Schensted (1960) that the set of such functions 
is complete for channel flows and consequently there is no conceptual diffi- 
culty in the prescription. Lastly, there is no continuous eigenspectrum for 
viscous channel flows. Thus, transient growth in these flows is due entirely 
to the non orthogonal eigenvectors of the non self-adjoint differential equa- 
tions. The optimal formulation is due to Farrell (1988), Butler and Farrell 
(1992) and has been corroborated by Reddy and Henningson (1993). Up 
to thirty eigenmodes were employed in these works with relative amplitudes 
ranging from order one to one thousand. Typical output was a display of the 
normalized kinetic energy of the perturbations. It is in terms of this quan- 
tity that shows explicitly that transient growth leads to the large amplitudes 
cited before there is eventual decay. The initial distributions used for the 
optimal values are given in Farrell (1988) and Butler and Farrell (1992). In 
terms of the initial- value problem this is important information and certainly 
bears heavily on the possible physical realization. Examination of the initial 
conditions that produce optimal growth show that the optimal initial con- 
ditions for strictly two-dimensional disturbances have discernible structure 
while those that produce optimal growth of three-dimensional disturbances 
are rather nondescript. Thus it should be relatively easy to choose initial 
conditions that approximate the three-dimensional optimal conditions while 
some care must be taken to choose initial conditions which will approxi- 
mate the optimal two-dimensional disturbance. This aspect of realization 
of the optimal conditions is a major point of this study. Butler and Farrell 
(1994) show that a threshold amplitude exists for optimally configured two- 
dimensional initial conditions. For amplitudes above the threshold, transient 
growth leads to a nonlinear evolution to quasi-steady finite amplitude struc- 
tures, and for amplitudes below the threshold, the decay rate for long time 
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behaviour is predicted by the slowest-decaying Orr-Sommerfeld mode. The 
optimization procedures using the eigenfunction expansions in these studies 
are dependent on the fact that the discrete eigenfunctions form a complete 
set and therefore are not directly applicable to a problem such as the bound- 
ary layer where the complete set was shown by Salwen and Grosch (1981) 
to be a linear combination of discrete and continuum eigenfunctions. Butler 
and Farrell (1992) did study the boundary layer optimization problem by 
using the channel flow solutions to represent the continuous spectrum as a 
discrete set of modes. The solution procedure chosen here does not present 
the same difficulties in generalization to a problem with a continuous spec- 
trum since it does not depend on the expansion of an initial condition in 
terms of its eigenfunctions. 

The presentation of Gustavsson (1991) was somewhat different. Here, 
it was considered sufficient to take only one of the normal modes from the 
Orr-Sommerfeld equation and combine this with six from the accompanying 
Squire equation. In effect, this is an initial-value problem but extremely 
limited. Similarly, a measure of the energy is displayed but with the concen- 
tration on the energy gained by the Squire mode. The unknown coefficients 
here were determined by requiring the Orr-Sommerfeld mode to have an en- 
ergy of unity at time zero. It was shown that, if there is sufficient obliquity 
of the waves (nearly 90 degrees), the energy for the Squire mode can tran- 
siently increase to almost a factor of 100. In this case, three-dimensionality 
is crucial but this in no way makes it accessible as a method for realizing 
large transient growth in an arbitrary initial-value problem. Gustavsson did 
consider the effects of symmetrical and asymmetrical interaction of the two 
sets of eigenfunctions. There are decidedly different responses depending 
upon which choice is made. In particular, a symmetric Orr-Sommerfeld 
mode interacting with an asymmetrical Squire mode leads to the largest 
amplitude for the energy. In every case, there is eventual decay after the 
maximum as time advances for subcritical values of the Reynolds numbers. 

Like the boundary layer, viscosity is critical to the perturbations in chan- 
nel flows. Without a viscous fluid, there is no instability in any of these 
flows w r hen cast in terms of normal modes. Case (1960) investigated plane 
Couette flow in an inviscid fluid and, by the use of Laplace transforms, 
demonstrated that there are no normal modes at all, whether decaying or 
growing. At the same time he also showed that there are other solutions 
to this problem if one makes the proper analysis. The additional solutions 
are those due to a continuous eigenspectrum. Although not covered to this 
extent, plane Poiseuille flow must follow in like manner. Once the fluid is 
viscous, then the results are reversed, i.e., there is no continuous spectrum 
for channel flows but it has been proven that there is an infinite number of 
normal modes and the set is complete (DiPrima and Habetler, 1969). The 
boundary layer (see Gustavsson 1979, e.g) possesses both spectra because it 
is a semi-infinite problem. Parenthetically, however, it is interesting to note 
that there are only a finite number of normal modes for the boundary layer 
(see Mack 1976; Grosch and Salwen 1978) with the number depending upon 
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the value of the Reynolds number. Salwen and Grosch (1981) showed how 
an arbitrary initial disturbance can be expanded in terms of the complete 
set of discrete and continuum eigenfunctions, but determining an optimal 
initial condition would be difficult. It is only for infinite Reynolds number 
that an infinite (all damped) set of normal modes is possible. 

Physically, viscosity causes instability in much the same way that a 
spring would be unstable if the spring force had a time delay. Mathe- 
matically, it is a question of phasing of the perturbations. This is best 
illustrated by considering the equation for the total kinetic energy of the 
perturbations. The time rate of change for this energy depends upon two 
essential terms: dissipation due to viscosity and production manifested by 
one of the Reynolds stress components interacting with the mean flow to 
transfer energy and over-come the viscous dissipation. In the case of chan- 
nel flows this process is possible for Poiseuille but not for Couette flow. And, 
the boundary layer is also unstable in these terms. 

Because of the Squire theorem and the fact that a stability boundary w r as 
thought to be sufficient, three-dimensionality per se was neglected because 
it was recognized that purely two-dimensional disturbances had the largest 
growth factors. Another often overlooked reason is the fact that the Orr- 
Sommerfeld equation is only fourth rather than sixth order. This result is 
fortuitous, makes for ease in the mathematics, but omits important physics. 
The link to resolution of this problem is the Squire mode equation that is 
coupled to the behavior of the Orr-Sommerfeld equation so long as there is 
three dimensionality. It is this pair of governing equations that has formed 
the bases for the cited studies in transient temporal behavior. Except for 
the investigation of Gustavsson (1979) that uses the Laplace transform for 
the boundary layer, the work has relied completely on solutions in terms of 
normal modes. If the transient dynamics is to be important, then the effects 
of specific initial conditions must be examined together with any optimum 
strategy. Certainly, three-dimensionality must be treated in a thorough 
manner. These aspects of the linear perturbation problem are the central 
goals of this presentation. This is done with the use of numerical integration 
of the governing partial differential equations. In no way is an expansion 
in normal modes suggested but, at the same time, it will be seen that all 
of the known results of classical stability theory as well as the optimization 
problem can be reproduced. It is further stressed how important are the 
details of the perturbation field, most notably the vorticity. 

2 BASIC GOVERNING EQUATIONS 


For plane viscous channel flows, the fluid is taken as incompressible with 
the basic flow parallel, U — U(y),V — W = 0. Then, the nondimensional 
linearized equations of motion can be written as 

du dv dw _ 

dx + dy dz ’ 
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( 2 ) 


du 

~di 


du dU dp _] 

+ v r, + ^ + ai = Re 


d 2 u d 2 u d 2 u 

dx 2 dy 2 + dz 2 ’ 


dv (h frp _ i (Pv_ (Pv d^v 
~di + 1 !h + dy dx 2 dy 2 dz 2 \ ’ 


(3) 


and 


dw TT dw dp _! d 2 w d 2 w d 2 w 
-r— + U -5 — |- -r— = Re 2 + ITT + A 2 • 

dt dx dz ydx 2 dy 2 dz 2 


(4) 


Here, Re = U 0 h/v is the Reynolds number, where h is the channel half- 
width, U 0 the centerline velocity and u the kinematic viscosity. Time is 
nondimensionalized by the advective time scale h/U 0 . On using the Fourier 
transformations defined with respect to x and 2 as 


v(a;y;r,t)= [ f v(x,y,z,t)e t{ax+lz) dxdz, 
J —OO J — OO 

equations (1) to (4) become 


and 


.. v dv 

— i{au + 7 w) + — = 0, 


du 

'dt 


— iaU u + U'v — iap = Re 1 


d 2 u 


dy 2 


— d 2 u 


dv dp , 

— iaU v + — — Re 

dt dy 


- iaU w — i'jp = Re 1 
dt 


<Pv 

dy 2 

d 2 w 


dy 2 


- 2 w 
- a v 


— d 2 w 


respectively, with U' = dU/dy and a 2 = a 2 + j 2 . 
The Squire transformation, written as 


(5) 

( 6 ) 

(7) 

( 8 ) 

(9) 


au + 7 w — du (10) 

— ju + aw = aw (11) 

and combined with operations on (6) to (9) enables us to obtain the pair of 
equations 




+ iaU"i> = Re 1 


d 4 v 2 d 2 # 


-I- d 4 v 


( 12 ) 


and 



w — sin pU'v + Re 1 


d 2 w 

dy 2 


— a 2 w 


(13) 


where sin <f> = 7/a and w is proportional to the normal vorticity component 
(u y = iaw). The first term on the right hand side of (13) is called the vortex 
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tilting term which acts as a forcing term to the normal vorticity component. 
The vortex tilting term is a product of the mean vorticity in the spanwise 
direction (cj z = - U * ) and the perturbation strain rate (dvfdz), and for 
a three dimensional disturbance, gives rise to the increase of the normal 
vorticity component. It is clear that the solutions of (12) and (13) combined 
with continuity and the Squire transformation are equivalent to solving (6)- 
(9). Likewise, p can be determined from (9). In either case, solutions of 
the equations are subject to imposed initial conditions and the following 
appropriate boundary conditions at the channel walls 

«(±M) = |j(±M) = «K±M) = 0- ( 14 ) 

For the mean velocity, we shall consider both plane Poiseuille flow 

U(y) = 1 - y 2 


and plane Couette flow 7 

U{y) = y- 

To evaluate the other velocity components, the quantities v and w are 
first computed from (12) and (13), respectively. Then the Squire transfor- 
mation (10)-(11) is inverted to give u and w 


v ? cos (bdv , . _ 

u - sin <pw, 

a ay 

v i sin <fidv 

w ~ — — + cos <pw . 

a ay 


(15) 

(16) 


By knowing the velocity components, the vorticity components can be deter- 
mined in a straightforward manner by appealing to their definitions, yielding 


Wt 


dw . „ 

a - + *7», 

oy 


Uy = —i^u + iaw = idni, 


(17) 

(18) 


and 


u, = — tav — 


du 

dy' 


(19) 


Finally, we remark here that, if one seeks solutions to (12) and (13) 
of the form then (12) becomes the more familiar Orr-Sommerfeld 

equation and (13) the Squire equation. Solutions of these equations will 
yield classical normal modes and normally (i) transient dynamics and (ii) 
effects of various initial conditions are ignored. At sub-critical Reynolds 
numbers, where the normal modes are damped, transient behavior may be 
extremely important. A variety of authors used eigenfunction expansions 
to examine transient dynamics (e.g., Gustavsson, 1991; Butler and Farrell, 
1992; Reddy and Henningson, 1993; as main references). In particular, 
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since the eigenfunctions form a complete set (DiPrima and Habetler, 1969; 
Herron, 1980), solutions to (12) and (13) were sought in the form 


w(y,t) 


= £ 


Aje-**** 


3 


vj(y) 

w ]{y) 


i 


0 

Wj {y) 


where {Aj} and {^j} are the eigenvalues of the Orr-Sommerfeld and Squire 
equations, respectively, with the eigenvalues distinct. As already noted, ei- 
ther (i) variational methods are then used to determine optimal growth; or 
(ii) a finite combination of eigenfunctions, though extremely limited, are 
then chosen and their subsequent transient behavior followed. This is not 
necessarily a weakness because an infinite set is available but only a finite 
number have been used. An alternative, yet novel, approach is to solve the 
system (12) and (13) directly by a simple numerical scheme. While this ap- 
proach does not directly select the optimal initial conditions that provide the 
optimal algebraic growth for a given set of parameters (d, (j>, Re), it does al- 
low one to follow rather easily the transient dynamics of any given prescribed 
initial condition and to determine if optimal growth can be approximated by 
realizable initial conditions. Also, simple maximization methods (see section 
4.3) can be easily applied to the numerical solution in order to select the 
optimal initial condition in a rational manner. Farrell and Moore (1992) also 
integrated the governing equations for oceanic flows, but again their focus 
was on determining an optimal initial condition by repeated integration of 
the perturbation equations and its adjoint and not on the dynamics of spe- 
cific initial conditions. Here we show that the optimal initial condition can 
be determined in a straight forward manner which circumvents the necessity 
of using the adjoint solution, is conceptionally easier to understand, and is 
easier to implement than the adjoint method or the eigenfunction expansion 
method. The approach, however, is not necessarily computationally faster. 
Furthermore, this approach is more robust than employing eigenfunction 
expansions since it can be applied to other classes of problems where only a 
finite number of normal modes exist; e.g., boundary layers, free shear layers 
(see Criminale, et al. 1994 for related work on inviscid flows). 


3 NUMERICAL SOLUTIONS 

The partial differential equations (12) and (13) were solved numerically 
by the method of lines. The spatial derivatives were center differenced on 
a uniform grid within the channel while one-sided differences were used at 
the walls. The resulting system was then integrated in time by a fourth- 
order Runge Kutta scheme with all calculations done in double precision. 
The results were checked for convergence by increasing the number of mesh 
points, varying between 500 mesh points for low Reynolds numbers to a 
maximum of 10,000 at larger Reynolds numbers. The table below shows the 
numerically computed growth rate for plane Poiseuille flow as a function of 
grid points for Re = 10,000, <5 = 1 and <f> = 0°: 


7 



GRID POINTS 

GROWTH RATE 

500 

0.003726 

1000 

0.003736 

2000 

0.003739 


The exact value from Orszag (1971) is 0.00373967. The number of grid 
points were sufficient to resolve the boundary layers near the walls. No ef- 
fort was made to optimize the number of grid points by employing nonuni- 
form meshes. If this were done, far less grid points would be needed. All 
calculations presented in this paper represent converged solutions. 

Before investigating the effects of various initial conditions and their 
subsequent transient behavior, it was first instructive to compare for plane 
Poiseuille flow numerically computed growth rates and eigenfunctions to 
those of the Orr-Sommerfeld at super-critical Reynolds numbers. Figure 
1(A) shows the growth rates obtained from the numerical solution of (12) 
(shown as circles) and those obtained from the Orr-Sommerfeld equation 
(shown as the solid curve) for Re = 10, 000 and <f> = 0°. The agreement is 
excellent. The corresponding real and imaginary parts of the eigenfunctions 
from both the numerical solution (solid) and the Orr-Sommerfeld solution 
(dashed) are displayed in figures 1(B) and 1(C), respectively, for the case 
a = 1 and <j> = 0°. Note that the two curves essentially lie on top of each 
other. To demonstrate the strength of the procedure, figure 2 shows similar 
results but at Re = 10 6 . 


4 PERTURBATION ENERGY 


As mentioned above, we are particularly interested in the effects of vari- 
ous initial conditions and their subsequent transient behavior at sub-critical 
Reynolds numbers. In order to examine the evolution of various initial 
conditions, the energy density in the (a,4>) plane as a function of time is 
computed. The energy density is defined as 


E(t; a, <f>. Re) = [ju 2 | + |u 2 | + |u; 2 |] dy. (20) 

The total energy of the perturbation can be found by integrating (20) over 
all a and <fi. A growth function can be defined in terms of the normalized 
energy density, namely 


C?(f ; d, <p, Re ) 


E{t\ a, 4 i, Re) 
E{ 0; d,<£, Re) 


( 21 ) 


measures the growth in energy at time t for a prescribed initial condition at 
t = 0 . 

Various initial conditions are used to explore transient behavior at sub- 
critical Reynolds numbers. For the integration of (12) the initial conditions 
for v are provided in the table below; Do — Dq(»,7). 
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CASE 

v(y, 0) 

I 

Oo(i - y 2 ) 2 

II 

^[cos (3 - cos ((3y)Y, [3 = mr 

III 

$[cos(3- cos((3y)] e J-^ ; (3 - mr 

IV 

toovO- - y 2 ) 2 

V 

ygf[cos/3 - cos(/fy)] e v/i ^ J ; (3 - mr 


Note that the first three cases correspond to symmetric initial conditions 
while the last two are asymmetric. For the integration of (13) the initial 
conditions for ii> are provided in the table below. 


CASE 

w(y, 0) 

i 

0 

ii 

fii cos (/3 0 y); Po = (2 n - l)i r/2 

iii 

Qi sin(/?iy); A = mr 


The second initial condition is symmetric while the last is asymmetric. We 
remark here that, in choosing these particular initial conditions, no attempt 
was made to find an optimal initial condition that correspond to a maximum 
transient growth. Also, particular eigenmodes are not investigated since 
these were studied previously by Gustavsson (1991). Rather, reasonable ini- 
tial conditions consisting of polynomials and transcendental functions were 
constructed and their subsequent transient behavior followed. The above 
set of initial conditions is general enough that it should be sufficient to test 
whether the optimal transient growth previously determined is realizable or 
not. 

4.1 TWO-DIMENSIONALITY 

In this section results are presented for plane Poiseuille flow with a = 
1.48, 0 = 0° and a Reynolds number of Re = 5,000, and for plane Couette 
flow with q = 1.21, 4> = 0° and Re = 1,000. For plane Poiseuille flow, 
this Reynolds number is sub-critical and thus the growth function G will 
eventually decay in time to zero. Because plane Couette flow is linearly 
stable, G will always eventually decay. These values correspond to those 
found by Butler and Farrell (1992) who showed that these values are the 
best 2-D optimal using variational techniques in the sense that they give the 
largest algebraic growth. 

The growth function is plotted as a function of time in figure 3 for plane 
Poiseuille flow and in figure 4 for plane Couette flow. In both figures, the 
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curves in (A) correspond to all possible combinations of the initial conditions 
for v and w with Q 0 = f} 1 = A = n = l. Except for some minor algebraic 
growth in the case of plane Couette flow, the growth function decays. This 
decrease in the amplitude of the maximum growth is somewhat unexpected 
since previous work shows that by considering optimal initial conditions that 
substantial growth can occur even for these two dimensional disturbances. 
Therefore, this issue was pursued further by considering initial disturbance 
velocity profiles with more zeroes which, in some sense, corresponds to the 
higher eigenmodes whose inclusion in the optimization analysis was neces- 
sary to achieving the high growth rates. Figures 3(B) and 4(B) are for the 
initial condition (II, i) w r ith various values of n. For plane Poiseuille flow and 
with n — 7 the maximum value is 12, and for plane Couette flow r and with 
n = 3 the maximum value is 4.8. In both cases, moderate transient growth 
is observed, with the maximum growth being lower than that obtained by 
Butler and Farrell (1992). It is interesting to note that an “optimal” initial 
disturbance will be found for both the Couette flow 7 and the Poiseuille flow 
that is consistent with the optimal conditions determined by Butler and Far- 
rell (section 4.3). For Couette flow, the maximum optimal energy growth for 
this choice of a and 4> occurs at t = 8.7. Here, we observe that the largest 
growth is for the initial condition with n = 3 and the maximum occurs at 
time t = 7.8. The same can be said of Poiseuille flow. The optimal initial 
conditions produce a maximum at time / = 14.1 and the largest growth here 
is for n = 7 that has a maximum at time t = 14.4. It is easy to see how r these 
solutions for different values of n can be combined to produce an optimal 
solution. This issue is explored further in section 4.3. 

4.2 THREE-DIMENSIONALITY 

In this section results are presented for plane Poiseuille flow 7 with d = 
2.044, <f> = 90° and a Reynolds number of Re = 5, 000, and for plane Couette 
flow with d = 1.66, <f) = 90° and Re = 1,000. These values appear in 
Butler and Farrell (1992) and the choice of a corresponds to the stream wise 
vortex with largest growth. For plane Poiseuille flow the global optimal 
coincides with the streamwise vortex but not so for Couette flow. In the 
latter case, the global optimal w r as shown to be at <fi « 88°. Here, we are 
only interested in presenting two representative cases. In all calculations we 
set flo = JIi = A = n = 1, Changing the sign of Sli produced only very 
minor changes in the solutions. Numerical experiments were also carried 
out for various combinations of n and A. In all cases, the largest value of 
G are for the cases presented below 7 , namely n = A = 1, and therefore only 
these cases are presented here. 

The growth function is plotted as a function of time in figure 5 for plane 
Poiseuille flow and in figure 6 for plane Couette flow. The curves in (A) cor- 
respond to the initial condition (i) for u), w r hile the curves in (B) correspond 
to (ii) and the curves in (C) correspond to (iii). Comparing figures 3 and 5 
for the case of plane Poiseuille flow 7 , and figures 4 and 6 for the case of plane 
Couette flow, we see that the transient growth is significantly larger for three 
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dimensional disturbances than those for two dimensional disturbances. For 
plane Poiseuille flow, the maxima of the symmetric disturbances (labeled 
I, II, and III) and the maxima of the asymmetric disturbances (labeled TV 
and V) in figure 5(A) are within 90% of the global maxima reported by 
Butler and Farrell (1992). They point out that the presence of streamwise 
vorticity, while passive to nonlinear dynamics (Gustavsson, 1991), can cause 
the development of streaks which may themselves be unstable to secondary 
instabilities or possibly produce transient growth of other types of pertur- 
bations. For plane Couette flow, the maxima of the symmetric disturbances 
(labeled I, II, and III) in figure 6(A) are within 97% of the maxima reported 
by Butler and Farrell, while the maxima of the asymmetric disturbances 
(labeled IV and V) are significantly smaller. The significance of this is that 
any initial condition with v velocity symmetric and no initial vorticity will 
give near optimum results when three-dimensionality is considered. This 
easily explains the growth observed by Gustavsson (1991) when a limited 
normal mode initial condition was chosen. 

Since this large transient growth for three dimensional disturbances is 
a direct result of the growth of vorticity, it is necessary to ask whether the 
growth is a special case because the initial vorticity is neglected, or will 
the growth in energy remain once the energy of non-zero initial vorticity 
is included in the calculation. Note in figures 5 and 6 that the symmetric 
disturbances all behave in a similar manner when the initial condition for 
w is given by (i), but differ substantially when conditions (ii) or (iii) are 
used. It is clear that the substantial growth in energies shown in figure 5(A) 
is directly attributable to the generation of normal vorticity through the 
coupling term in equation (13). This produces growth factors of 3900-4200 
for symmetric initial profiles (I, II, and III) and growth factors of 2200-2400 
for anti-symmetric initial profiles (IV and V). The responses change signif- 
icantly when the energy of the initial profiles (ii) and (iii) are included in 
the initial normalization of the growth factor G. The responses to initial 
profiles I and IV are lowered but still show significant growth. The initial 
energy of II and IV are larger than those of (ii) and (iii). However, when the 
energies of the initial normal vorticities (ii) and (iii) are significantly larger 
or comparable to the energies of the initial velocity profiles II, III, and V, the 
substantial transient growth as measured by the total energy can decrease 
by a factor of ten. This is not to say that the transient growth shown in (A) 
is not also present in (B) and (C), but that the normalization in (B) and (C) 
reflects a more proper measure of growth that includes the total energy of 
an arbitrarily chosen initial disturbance and not just an optimal one. This 
suggests that the algebraic growth as measured here and in previous work 
is extremely sensitive to the presence of any initial normal vorticity, specif- 
ically to the inclusion of the energy associated with the initial streamwise 
vorticity in normalization of the growth factor. Vorticity components are 
presented in figure 7 for plane Poiseuille flow with initial conditions given by 
(II, i) and (II, ii), respectively, showing that the maxima in time of the nor- 
mal and spanwise vorticity components decrease when the normal vorticity 
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component is initially non-zero. The streamwise component decays and is 
unaffected by this change of initial conditions. 


Table 1: Maximum values at various Re corresponding to the initial condi- 
tion (I,i) for plane Poiseuille flow with a = 2.044 and 4> = 90°. 


Re 

t* . 

G* 

t* j Re 

G*/ Re 2 

500 

41.0 

45.6 

0.0820 

1.824xl0 -4 

1000 

82.0 

181.9 

0.0820 

1.819xl0 -4 

2000 

164.0 

726.9 

0.0820 

1.817xl0 -4 

4000 

327.0 

2907.0 

0.0818 

1.817xl0 -4 

5000 

409.0 

4542.0 

0.0818 

1.817xl0 -4 


Table 2: Maximum values at various Re corresponding to the initial condi- 
tion (I,i) for plane Couette flow with a = 1.66 and <j> = 90°. 


Re 

t* 

G* 

R/Re 

G*/Re' 2 

500 

69.0 

291.4 

0.1380 

1.165xl0 -3 

1000 

139.0 

1165.2 

0.1390 

1.165xl0 -3 

2000 

277.0 

4660.7 

0.1385 

1.165xl0~ 3 

4000 

554.0 

18642.4 

0.1385 

1.165x10^ 

8000 

1109.0 

74569.6 

0.1386 

1.165xl0 -3 


Finally, shown in Tables 1 and 2 are the times (£*) for which the growth 
function attains its maximum (G*) at various Reynolds numbers for plane 
Poiseuille and plane Couette flow, respectively. Note that the time scales 
as Re while the maximum scales as Re 2 , as previously pointed out by Gus- 
tavsson (1991). 

4.3 OPTIMAL INITIAL CONDITIONS 

A mechanism for rapid transient growth when the initial conditions is 
expressed as a sum of the eigenfunctions has been explained by Reddy and 
Henningson (1993). The concept is that a group of eigenfunctions are nearly 
linearly dependent so that, in order to represent an arbitrary disturbance 
(say), then it is possible that the coefficients can be quite large. Now, since 
each one of these nearly linearly dependent eigenfunctions have differing 
decay rates, the exact cancellations that produce the given initial distur- 
bance might not persist in time and thus significant transient growth can 
occur. This process can (and is) taken a step further in order to determine 
the optimal initial condition (still expressed as a sum of the non-orthogonal 
eigenfunctions) that produces the largest relative energy growth for a certain 
time period. This process is completed and it does have the feature that the 
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nearly linearly dependent eigenfunctions are multiplied by coefficients three 
orders of magnitudes greater than the others. This optimal initial condition 
produces a growth factor of about 20 for the two-dimensional disturbance 
in Poiseuille flow. However, this optimal growth is nearly destroyed by not 
including the first eigenfunction (growth drops to a factor of 6 rather than 
twenty) which seems to indicate that the prior explanation of (initial) ex- 
act cancellations by the nearly linearly dependent eigenfunctions is not the 
entire mechanism. Butler and Farrell (1992) also calculated optimal initial 
conditions in terms of a summation of the eigenfunctions (although they 
put no particular emphasis on the importance of using this approach) and 
re-iterated the importance of nearly linear dependence of the modes to the 
transient growth. Butler and Farrell (1992) also explained the transient 
growth of the optimal initial conditions in terms of the vortex-tilting mech- 
anism and the Reynolds stress mechanism, since these (physical) arguments 
apply no matter what the solution method; perhaps in the end they should 
be preferred to those based strictly on the mathematical procedure, espe- 
cially since owing to the non-orthogonality and near linear dependence of 
the eigenfunctions. The resulting optimal initial disturbance is difficult to 
visualize physically for someone who is not privy to the calculations. 

Using the method that we have been following, an optimization proce- 
dure can be determined without resorting to a variational procedure requir- 
ing the numerical determination of the eigenfunctions. Now for these channel 
flows, the difficulties of calculating the eigenvalues and eigenfunctions are 
well known. A group of eigenfunctions are nearly linearly dependent and the 
whole set is not orthogonal. Also, the issue of degeneracies in the parameter 
space must be considered. Thus, an attempt to write a specific initial con- 
dition as an eigenfunction expansion (which also requires the determination 
of the eigenfunctions to the adjoint equation) runs into difficulties for the 
uninitiated. In this paper, we show that calculating growth for particular 
initial conditions does not present great difficulty. Furthermore, a closer 
inspection of initial condition II suggests that each of these disturbances is 
in essence a single Fourier mode of an arbitrary initial condition. If one 
were to consider an arbitrary odd function for the u velocity satisfying the 
boundary conditions written in terms of a Fourier sine series, then the initial 
condition in the v velocity is given by II. Thus, if one wished to determine 
an optimal initial disturbance, a maximization procedure could be applied 
to an arbitrary linear combination of these modes, all of which are initially 
orthogonal and linearly independent. Clearly, if one wanted to also include 
non-zero initial vorticity in such an optimization scheme it would not be 
difficult to include (and these initial conditions are of course very important 
when modeling real disturbances as opposed to optimal disturbances). The 
results presented here show that, if included in the optimization procedure, 
the initial vorticity modes would not contribute to the optimal solution for 
the cases considered. 
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To start our optimization scheme, we consider the total solution u = (u, v, w) 
to be the sum 


N 


<) = + ib k)uk{y,t) 

k=i 


where each of the vectors u k (y, t) represents a solution to equations (12) and 
(13) subject to the initial conditions 


u(y, 0) 


ia 
ki r 


cos cf) sin k'Ky 
(cos kit — cos kity) 
sin 0 sin kity 


In order to maximize the growth function, it is sufficient to maximize the 
energy ^ 

E(t) = J ^u(y,t)-u*(y,t)dy 
subject to the constraint 

E( 0) = 1. 


Therefore, we use Lagrange multipliers to maximize the function 


<?(<) = J u{y,t)-u*{y,t)dy- X (^ ^ u(y,0) • u*(y,Q)dy - l) 


which requires 


dG 

0a k 


= 0 , 


db k U) 


k = 1,2, ...,1V. 


The set of equations thus derived produce a 2 N x 2 N generalized eigenvalue 
problem. A search over the eigenvectors gives the initial condition with 
initial unit energy that maximizes the function G at time t. 

To illustrate the optimization procedure, we perform the calculations for 
the two cases reported by Butler and Farrell (1992). The first is the compu- 
tation of the two-dimensional optimal for a = 1.48, <f> = 0° and Re = 5,000. 
In figure 8(A) we show the growth factor at t = 14.1 for each individual mode 
as well as for the optimal solution for various values of A T . The convergence 
as N oc is nicely illustrated; compare to Reddy and Henningson (1993), 
e.g. In figure 8(B) the magnitudes of the coefficients which produce the op- 
timum with N — 20 are shown. There are no surprises. Each coefficient is of 
reasonable size, with the largest coefficient being a factor of ten greater than 
the first coefficient, and not a factor of a thousand as is the case when using 
eigenfunction expansions. The magnitudes peak for n = k, 7 and 8 which 
could be easily predicted from the previous graphs for the responses to each 
individual mode. For completeness the initial velocity profiles are shown in 
figure 9. These are consistent with the the initial perturbation streamfunc- 
tion contour plots shown in Butler and Farrell (1992). A similar calculation 
could be made for Couette flow but is unnecessary since the relevant infor- 
mation is easily determined by examination of figure 4(B). The magnitudes 
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of the coefficients which produce the optimal growth at t = 8.7 peaks be- 
tween modes 3 and 4, and converge quickly as N -> oc. The resulting initial 
condition is consistent with the initial perturbation streamfunction contour 
plots shown in Butler and Farrell (1992) for the same case. The second cal- 
culation is for the optimal three-dimensional disturbance. The parameters 
chosen are Re = 5,000, d = 2.044 and <j> = 90°. The initial conditions that 
produce a maximum growth at t = 379 are found. The results are shown 
in figure 10, and the composition of the initial conditions in terms of the 
modes chosen here could be easily determined from the individual responses 
of each mode. The key to this observation is that the chosen modes are 
not nearly linearly dependent as are the eigenfunctions and indeed provide 
a rational and easily understood basis for the calculation of arbitrary initial 
conditions. 

It must be re-iterated that, although it is possible and conceptually easy 
to reproduce the optimal initial conditions that have been previously found, 
the maximum transient growth is only a measure of what is possible and 
not what will actually occur as has been the difficulty in experiments. It is 
at least as important to investigate whether such large growth is possible 
for arbitrary initial conditions. In this regard, the results presented here 
produce a mostly negative answer to this question. For two-dimensional 
disturbances in Poiseuille flow, the transient growth observed for arbitrarily 
chosen initial conditions using this approach are, at best, only 25% of the 
optimal. When considering a fixed wavelength a and a fixed obliqueness <j > , 
it is seen that very large relative energy growth of the perturbation can be 
observed in Poiseuille flow for oblique disturbances with arbitrary velocity 
profiles restricted to having zero initial normal vorticity, but the relative 
energy growth quickly decreases when arbitrary disturbances are combined 
with initial normal vorticity. Similar results are found for Couette flow. 

5 CONCLUSIONS 

Plane Poiseuille flow and plane Couette flow in an incompressible fluid 
have been investigated subject to the influence of small perturbations. In- 
stead of using the techniques of classical stability analysis or the more re- 
cent techniques involving eigenfunction expansions, the approach has been 
to first Fourier transform the governing disturbance equations in the stream- 
wise and spanwise directions only and then solve the resulting partial dif- 
ferential equations numerically by the method of lines. Unlike traditional 
methods where traveling wave normal modes are assumed for solution, this 
approach offers a means whereby completely arbitrary initial input can be 
specified without having to resort to the complexity of eigenfunction ex- 
pansions. Thus, arbitrary initial conditions can be imposed and the full 
temporal behavior, including both early time transients and the long time 
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asymptotics, can be determined. All of the stability data that are known for 
such flows can be reproduced. Finally, an optimization scheme is presented 
using the orthogonal Fourier series and all previous results using variational 
techniques and eigenfunction expansions are reproduced. 

The benefit of this novel approach is clear: it can be applied to other 
classes of problems where only a finite number of normal modes exist, such 
as the boundary layer. In addition, this numerical approach has recently 
been successfully applied to free shear flows in an inviscid fluid (Criminale 
et al., 1994). These concepts are being extended to the Blasius boundary 
layer in an incompressible or compressible medium. 
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(B) 



(C) 


Figure 1: (A) Plot of growth rates as a function of wavenumber a. The 
circles corresponds to the numerically computed values from the partial 
differential equation, and the solid curve corresponds to the growth rate 
computed using the Orr-Sommerfeld equation. (B) The real part and (C) 
the imaginary part of the eigenfunction as a function of y for a — 1. Results 
for plane Poiseuille flow with <f> = 0° and Re = 10 4 . 
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Figure 2: (A) Plot of growth rates as a function of wavenumber a. The 
circles corresponds to the numerically computed values from the partial 
differential equation, and the solid curve corresponds to the growth rate 
computed using the Orr-Sommerfeld equation. (B) The real part and (C) 
the imaginary part of the eigenfunction as a function of y for a = 0.44. 
Results for plane PoiseuiUe flow with <j> = 0° and Re = 10 6 . 


19 





Time 


Figure 3: Plot of the growth function G as a function of time for various 
initial conditions. Results for plane Poiseuille flow with a — 1.48, 4> = 0°, 
and Re — 5, 000. 
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Figure 5: Plot of the growth function G as a function of time for various 
initial conditions. Results for plane Poiseuille flow with a = 2.044, 4> = 90°, 
and Re = 5, 000. 
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Figure 6: Plot of the growth function G as a function of time for various 
initial conditions. Results for plane Couette flow with a = 1.66, <f> ~ 90°, 
and Re — 1,000. 
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(A) 



(B) 


Figure 7: Plot of the vorticity components as a function of time for the 
initial conditions (A) (II, i) and (B) (II, ii). Results for plane Poiseuille flow 
with a = 2.044, d> = 90°, and Re = 5,000. 
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Figure 8: (A) Plot of the growth function G at t = 14.1; individual mode 
results denoted by *; cumulative results from optimzation procedure denoted 
by o. (B) Plot of the magnitude of the coefficients ak + ibk from optimzation 
procedure for N = 20. For plane Poiseuille flow with a = 1.48, <f> — 0°, and 
Re = 5,000. 
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Figure 9: Plot of velocities for plane Poiseuille flow from optimzation pro- 
cedure with N = 20, a = 1.48, <j> = 0°, and Re = 5,000. (A) u T (solid), 
Ui (dashed) at t = 0; (B) v r (solid), V{ (dashed) at t = 0; (C) u r (solid), U{ 
(dashed) at t = 14.1; (D) v r (solid), Vi (dashed) at t = 14.1. 
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Figure 10: (A) Plot of the growth function G at t = 379.0; individual mode 
results denoted by *; cumulative results from optimzation procedure denoted 
by o. (B) Plot of the magnitude of the coefficients a * + ibk from optimzation 
procedure for N — 20. For plane Poiseuille flow with ct = 2.044, <f> = 90°, 
and Re = 5, 000. 
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